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ABSTRACT 


¥ 


A  combined  direct/ inverse  three-dimensional 
transonic  wing  design  method  is  presented.  The  method 
is  built  around  the  ZEBRA  II  transonic  potential  flow 
solution  algorithm  to  provide  a  design  method  that  is 
particularly  suited  for  use  on  a  vector  computer.  The 
development  of  a  pilot  design  computer  code  and  a 
baseline  design/ analysis  code  is  described.  Results 
are  presented  that  verify  the  accuracy  and  consistency 
of  the  design  method. 


ADMINISTRATIVE  INFORMATION 


The  work  presented  was  a  joint  effort  by  Lockheed-Georgia  Company  and 
Texas  A&M  University  supported  by  the  Naval  Air  Systems  Command  under  the 
cognizance  of  D.  G.  Kirkpatrick  (NAV AIR-3 1 ID ) ,  Navy  Contract 
N00167-81-C-0078-P00001.  The  authors  acknowledge  Dr.  Tsze  C.  Tai,  contract 
monitor  at  the  David  Taylor  Naval  Ship  Research  and  Development  Center,  and 
Mr.  Jerry  South,  NASA-Langley  Research  Center,  for  providing  computer  time. 


INTRODUCTION 

In  recent  years,  the  increasing  importance  of  transonic  flight  by  both 
military  and  commercial  aircraft  has  prompted  a  large  amount  of  research  to 
develop  more  accurate  and  reliable  computational  methods  for  the  analysis 
of  aircraft  configurations  in  transonic  flow.  This  research  was  spurred  by 
the  Increasing  costs  of  wind  tunnel  tests  and  the  Interference  and  scale 
problems  associated  with  tests  conducted  at  transonic  conditions.  As  a 
result  of  this  effort,  several  computer  codes  have  been  developed  to  calcu¬ 
late  the  transonic  flow  about  wing  and  wing- body  configurations.  A  few  of 
these  codes  have  demonstrated  levels  of  accuracy  and  reliability  that  have 
gained  them  acceptance  in  the  aircraft  industry  as  useful  analysis  tools. 

Unfortunately,  the  development  of  wing  design  codes  has  lagged  the 
development  of  analysis  codes.  The  research  undertaken  in  the  past  few 
years  to  develop  more  efficient  wing  design  methods  has  centered  on  two 
different  approaches  to  the  design  problem— numerical  optimization  and 
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inverse  design.  Numerical  optimization  provides  a  means  of  automating  the 
trial-and-correction  design  process  using  analysis  methods.  In  theory, 
optimization  allows  the  designer  to  specify  a  quantity  to  be  minimized, 
such  as  drag,  without  prior  knowledge  of  the  flow  details  that  will  produce 
the  objective  design.  However,  most  (if  not  all)  transonic  codes  cannot 
predict  drag  accurately  enough  to  use  it  as  a  design  objective  in  the  opti¬ 
mization  process.  Therefore,  the  difference  between  a  specified  design 
pressure  distribution  and  the  computed  pressure  distribution  at  a  span 
station  is  used  as  the  function  to  be  minimized  in  the  optimization  pro¬ 
cess.  However,  the  amount  of  computer  time  required  for  optimization 
limits  the  technique  to  performing  the  design  at  one  span  station  at  a 
time.  In  addition,  a  considerable  amount  of  computer  expertise  is  required 
to  effectively  implement  the  optimization  procedure. 

In  the  inverse  approach,  the  wing  geometry  is  computed  by  specifying  a 
desired  pressure  distribution  over  a  part  of  the  wing  and  then  solving  a 
mixed  Neumann  and  Dirichlet  boundary  value  problem  by  finite  difference 
techniques.  Since  more  than  one  span  station  at  a  time  can  be  designed  by 
the  Inverse  technique,  it  would  appear  to  be  much  simpler  to  use  and  cost 
less  than  the  numerical  optimization  procedure. 

The  present  combined  direct/ inverse  transonic  wing  design  program  is  a 
joint  effort  of  the  Lockheed-Georgia  Company  and  Texas  A&M  University 
(TAMU)  to  develop  an  inverse  wing  design  method  that  Incorporates  the 
latest  advances  in  computational  transonic  aerodynamics.  The  research 
program  includes  three  major  tasks: 

1.  Formulation  of  the  Design  Method. 

2.  Development  of  a  Three-Dimensional  Inverse  Pilot  Code. 

3.  Development  of  a  Baseline  Unified  Design/ Analysis  Code. 

The  formulation  of  the  inverse  design  scheme  and  the  development  of  a 
pilot  code  to  validate  the  design  method  was  conducted  at  Texas  A&M 
University.  This  pilot  code  is  based  on  the  ZEBRA  II  three-dimensional 
transonic  potential  flow  code  developed  at  NASA  Langley  by  South  et  al.1'^ 
The  development  of  the  baseline  unified  design/ analysis  code  was  performed 
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by  Lockheed-Georgla  Company  in  two  stages.  A  three-dimensional  analysis 
code  capable  of  solving  the  flow  about  swept,  tapered  wings  without  twist 
was  developed  first.  This  code  served  as  the  baseline  code  for  the 
development  of  the  unified  design/ analysis  code.  The  inverse  design  method 
developed  at  TAMU  was  then  implemented  into  this  analysis  code. 


FORMULATION  OF  THE  DESIGN  METHOD 

The  principal  goal  of  this  research  program  was  to  develop  a  design 
method  that  is  accurate,  fast,  and  economical  to  use.  To  meet  this  goal, 
it  was  decided  that  the  design  method  should  incorporate  the  following 
features: 

1.  The  inverse  scheme  would  be  based  on  the  direct/ inverse 
approach  developed  by  Carlson  at  TAMU  for  a  transonic 
airfoil  design. 

2.  The  potential  flow  solver  would  use  the  conservative 
form  of  the  full  potential  equation. 

3.  Wing  surface  boundary  conditions  would  be  applied  on  a 
mean  plane  in  a  Cartesian  grid  system  to  simplify  grid 
generation  and  wing  shape  calculation. 

4.  A  fast,  vectorlzable  solution  algorithm  such  as  approxi¬ 
mate  factorization  or  the  ZEBRA  II  algorithm  would  be 
used  in  the  potential  flow  solver. 

In  the  direct/ inverse  design  method,  the  leading  edge  geometry  of  the 
airfoil  (usually  the  forward  10  percent)  is  specified;  the  remaining 
portion  of  the  airfoil  is  computed  for  a  specified  pressure  distribution. 
This  eliminates  the  need  to  specify  a  boundary  condition  in  the  leading 
edge  stagnation  region. 

The  accuracy  of  transonic  flow  solutions  for  arbitrary  swept  wings 
depends  on  the  form  of  the  governing  equation  (i.e.,  full  potential  or 
small  disturbance),  the  finite  difference  scheme,  and  the  computational 
mesh  system  employed  in  the  solution  algorithm.  The  conservative  form  of 
the  full  potential  equation  provides  the  most  accurate  solution  for  highly 
swept  wings.  In  addition,  solution  of  the  conservative  full  potential 
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equation  ensures  that  the  condition  of  aero  mass  flux  across  surface 
streamlines  will  be  satisfied  for  inverse  design  cases.  This  condition  is 
critical  for  the  accurate  calculation  of  wing  shape. 


The  selection  of  a  Cartesian  computational  grid  system  in  lieu  of  a 
body-fitted  system  such  as  those  used  by  Jameson^  and  Holst^~^  was  based 
on  results  obtained  by  Purcell  and  Carlson^  for  two-dimensional  transonic 
flow.  Purcell  and  Carlson^  showed  that  sufficient  accuracy  can  be  obtained 
for  full  potential  equation  solutions  by  applying  the  full  surface  boundary 
condition  on  a  mean  plane  in  a  Cartesian  grid  system.  This  plane  can  be 
located  on  a  grid  line  or  situated  between  two  adjacent  grid  lines.  Two 
problems  are  avoided  by  using  the  Cartesian  grid  system  and  mean  plane 
boundary  conditions  for  inverse  design  calculations.  First,  the  computa¬ 
tional  grid  does  not  have  to  be  recomputed  each  time  the  wing  shape  is 
computed  in  an  inverse  design  case.  Second,  intermediate  calculation  of 
wing  shape  during  the  potential  flow  solution  is  avoided.  The  new  wing 
shape  is  computed  only  after  the  potential  flow  solution  has  converged  to  a 
desired  value. 

In  order  for  any  inverse  scheme  to  be  cost  effective,  the  potential 
flow  solver  must  be  fast  and  reliable.  In  addition,  the  solution  algorithm 
should  be  amenable  to  vectorization  for  use  on  current  supercomputers  such 
as  the  CYBER  205  and  CRAY  I.  Two  existing  algorithms  meet  these  require¬ 
ments:  the  AF2  scheme  developed  by  Holst4-^  and  the  ZEBRA  II  algorithm 
developed  by  South  et  al.  After  an  unsuccessful  attempt  to  Implement 
the  AF2  algorithm  using  the  Cartesian  grid  system  described,  the  ZEBRA  II 
scheme  was  selected  for  the  potential  flow  solver.  The  selection  of  the 
ZEBRA  II  scheme  proved  fortuitous  because  it  allowed  the  use  of  the  pilot 
code  developed  at  NASA  Langley  as  a  test  bed  for  developing  the  inverse 
design  scheme.  Since  this  code  also  used  a  Cartesian  mesh  system,  the 
Inverse  schemes  developed  at  TAMU  using  the  ZEBRA  II  code  could  be 
Implemented  directly  into  the  baseline  analysis/design  code  being  developed 
at  Lockheed. 


DEVELOPMENT  OF  THE  THREE-DIMENSIONAL  INVERSE  PILOT  CODE 


The  development  of  the  inverse  design  method  occurred  in  two  phases. 
In  the  first  phase,  a  scheme  based  on  the  small  disturbance  approximation 
to  the  surface  boundary  condition  was  developed.  In  the  second  phase,  this 
technique  was  extended  to  use  the  full  surface  boundary  condition  applied 
on  a  mean  plane  in  the  Cartesian  grid  system.  The  small  disturbance  code 
was  developed  first  because  the  NASA  ZEBRA  II  code  used  the  small  disturb¬ 
ance  approximation  of  the  surface  boundary  condition.  In  addition,  it  was 
felt  that  the  schemes  developed  for  the  small  disturbance  boundary  condi¬ 
tion  provided  a  logical  foundation  for  building  the  full  potential  scheme. 
Because  of  its  importance  in  the  development  of  the  inverse  design  scheme, 
the  ZEBRA  II  algorithm  is  described  in  detail  in  Appendix  A.  A  complete 
description  of  the  development  of  the  inverse  method  is  given  in  Reference 
8. 


SMALL  DISTURBANCE  DESIGN  METHOD 


The  characteristics  of  the  inverse  design  method  in  a  Cartesian  grid 
system  are  affected  by  the  placement  of  the  ZK)  plane  on  which  the  surface 
boundary  conditions  are  applied.  In  the  Langley  ZEBRA  II  code  the  Z-0 
plane  is  located  between  two  grid  lines.  The  surface  boundary  condition  is 
implemented  by  replacing  the  difference  approximation  for  4>z  in  Eq.  (A-3) 
of  Appendix  A  on  the  plane  KWNGT-1/2  or  KWNGB+1/2  as  shown  in  Figure  1, 
with  the  small  disturbance  approximation 


-  u 

•  dx 


(1) 


In  this  way,  the  surface  boundary  conditions  can  be  implemented  in  the 
solution  algorithm  without  using  dummy  values  of  potential  or  costly  inter¬ 
polations  from  the  actual  body  surface.  In  addition,  the  complexity  of  the 
computer  program  is  reduced. 


In  the  inverse  design  method,  the  Neumann  surface  boundary  condition 
is  replaced  by  a  Dirichlet  condition  in  which  the  potential  is  specified 
directly  as  a  function  of  a  desired  pressure  distribution.  In  the  present 
research,  procedures  were  developed  for  implementing  the  inverse  boundary 
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conditions  that  would  not  destabilize  the  convergence  of  the  potential  flow 
solution  and  would,  at  the  same  time,  maintain  the  features  of  the  ZEBRA 
algorithm  that  make  it  vectorlzable.  The  following  scheme  was  developed 
for  the  Cartesian  grid  system  used  in  this  research. 

Referring  to  Figure  1,  the  small  disturbance  approximation  of  the 
pressure  coefficient,  Cp  “  -2$x,  at  the  mid-cell  point  X  on  the  mean  plane, 
Z“0,  can  be  written 


°P  *  "  AX  **R  “  *1^ 

The  objective  is  to  compute  the  value  of  $  at  point  A  as  a  function  of  the 
pressure  specified  at  point  X.  This  is  accomplished  by  first  computing  <j)R 
and  ^  by  extrapolation  from  the  points  above  the  mean  plane.  Three  point 
extrapolation  yielded  the  most  accurate  results.  Therefore,  can  be 
written  as 


<t>R  *  Ma  +  +  C<f>c 


(3) 


where  for  an  evenly  space  grid 

A 


Vc 


(2A  * 


ZB)  (zA  ' 


ZCJ 


B 


z  z 
AZC 


(ZB  " 


ZA>  (ZB  ‘ 


ZC} 


(4) 


(ZC  ~  ZAJ  (ZC  ~  ZB) 

Substituting  Eq.(3)  and  a  similar  one  for  into  Eq.  (2),  and  solving  for 
$A  yields 


*A  "  *D  '  A  [B(*B  '  V  +  C(*C  *  V  +  T  \]  (5) 

When  sweeping  through  the  grid  in  the  streamwlse  direction,  the  poten¬ 
tial  at  the  points  corresponding  to  Point  A  at  each  inverse  station  in  the 
cross-plane  is  determined  using  the  previously  described  approach.  Notice 
that  this  approach  uses  "old"  values  at  the  i  cross-plane  and  "new"  values 
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at  1-1.  The  cross-plans  Is  than  solved  using  the  standard  ZEBRA  approach. 
Since  the  ZEBRA  scheme  solves  for  only.  It  does  not  know  which  points 
are  inverse  and  which  are  direct.  At  the  end  of  the  double  pass  ZEBRA 
loop,  all  points  are  updated  by  A$  including  the  points  corresponding  to 
point  A.  To  correct  point  A  an  additional  calculation  is  performed  to  get 
back  to  its  boundary  value,  l.e.  ■  ($A  +A4>)-A$.  This  approach  is 
needed  to  retain  the  vectorization  feature  and  associated  efficiency  of  the 
ZEBRA  scheme.  The  actual  ZEBRA  loop  remains  blind  to  whether  the  station 
is  inverse  or  direct.  The  boundary  condition  alone  is  changed. 


After  a  converged  inverse  solution  is  obtained,  the  airfoil  shape  can 
be  determined  by  integrating  the  airfoil  slopes  obtained  from  the  wing 
boundary  condition,  l.e. 


dz 

dx 


(6) 


where  $z  must  be  obtained  from  the  inverse  solution.  The  first  attempts  to 
compute  $  at  the  wing  slit  used  the  4>  values  from  the  inverse  solution  and 
the  three  point  extrapolation  formulas.  However,  this  procedure  did  not 
yield  accurate  slopes  and  led  to  erroneous  airfoil  ordinates. 

A  second  approach  was  devised  that  used  the  finite  difference  approxi¬ 
mation  of  the  full  potential  equation  to  obtain  the  wing  slopes.  Expanding 
Eq.  (A-3 )  at  point  A  and  solving  for  $z  at  R  yields 

% '  £  ^ +  v»y  n  <7> 

R  i-A+4  y  J-A+*s 

The  value  of  wing  slope  at  each  inverse  station  is  then  obtained  by 
substituting  Eq.  (7)  for  $z  in  Eq.  (6).  Straightforward  trapezoidal 
integration  is  then  used  to  obtain  the  airfoil  ordinates  at  each  inverse 
station,  i.e.. 


i+1 


Ax  r,dz 


+  ir  ISr)  + 


(8) 
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FULL  POTENTIAL  DESIGN  METHOD 


As  In  the  small  disturbance  method,  the  full  potential  design  method 
uses  a  specified  pressure  coefficient  to  define  an  inverse  boundary 
condition.  However,  the  full  potential  pressure  coefficient  equation  is 
used  in  place  of  the  linearized  small  disturbance  equation.  The  full 
potential  pressure  coefficient  can  be  written  as 


C 


P 


— —  [  1  +  ^  ~  ^  N2 

,M2  2 


+  W* 


1] 


(9a) 


where  is  the  freestream  Mach  number,  Y  is  the  ratio  of  specific  heats,  q„ 
is  the  magnitude  of  the  freestream  velocity  vector,  and  u,v,w  are  the  local 
components  of  velocity  given  by  the  expressions 


u 


1  +  $ 

x 


Solving  Eq.  (9)  for  yields 


(9b) 


1  - 


(y  -  1)M2 


[('•  4s ) 


Hi 

Y 


-  1 


1  +  (?  +  o 


-  1 


(10) 


Referring  to  Figure  1,  this  expression  Is  applied  at  the  point  X  to 
extract  a  value  for  the  potential  at  point  A  as  was  done  in  the  small  dis¬ 
turbance  method.  However,  values  of  $ and  4>z  must  now  be  calculated.  As 
in  the  small  disturbance  method  a  three  point  Lagrangian  extrapolation  is 
used  to  define  values  of  $  on  the  wing  mean  plane. 

The  term  is  given  in  the  half  plane  at  Point  X  by  the  average  of 
at  points  R  and  L,  i.e.. 


1  ( JL_ 

2  v  2  Ay 


M  ■ 


The  value  of  4Z  at  point  X  is  computed  by  averaging  <f>z  at  L  and  R.  The 
expression  for  $z  in  the  half  plane  is  given  by  differentiating  the  general 
three-point  Lagrangian  extrapolation  with  respect  to  Z  so  that  4>z  at  R  is 
given  by 


t>2  -  A<J>a  +  B<j>B  +  C$c 


where 


-zB  ~2C 

<ZA  *  ZB>  (zA  -  ZC> 


— z  — z 
A  C 


<ZB  "  ZA>  (ZB  ’  ZC> 


~ZA  ~ZB 

(ZC  "  ZA}  (ZC  ’  V 


It  can  be  seen,  however,  that  evaluating  $z  and  $y  with  these  express¬ 
ions  poses  several  problems  when  used  in  Eq.  (10).  First,  Eq.  (10)  cannot 
be  evaluated  explicitly  at  each  time  step  since  u  appears  in  the 
denominator  of  the  right-hand  side  of  the  equation.  Second,  $a  appears  in 
both  and  expressions.  These  problems  were  overcome  by  updating  v/u 
and  w/u  only  every  ten  iterations  using  the  current  value  of  $A.  This  pro¬ 
cedure  is  based  on  the  fact  that  $y  and  $z  should  be  on  the  same  order  of 
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magnitude  as  the  slopes  and,  therefore,  very  small  relative  to  Thus 

the  denominator  in  Eq.  (10)  should  remain  on  the  order  of  one  in  the 
inverse  region. 


As  in  the  small  disturbance  method,  the  final  expression  for  $A  is 
given  by 


A 


1  - 


(Y  ~  DM2 


1  -  <J)2  +  (*)2 


Ax 

A 


(14) 


-i  (♦ 

A 


0  “  A 


where  A,  B,  and  C  are  the  Lagrangian  coefficients  given  in  Eq.  (4).  This 
boundary  condition  is  implemented  in  the  ZEBRA  algorithm  in  the  same  manner 
as  the  small  disturbance  boundary  condition.  The  ZEBRA  code  was  also  modi- 
fled  to  use  the  full  surface  boundary  condition  applied  on  the  wing  mean 
plane  for  calculations  in  the  direct  region. 

Calculation  of  wing  slopes  in  the  full  potential  method  is  also  per¬ 
formed  In  the  same  manner  as  was  done  in  the  small  disturbance  method  by 
using  the  residual  equation  to  define  the  w  velocity  on  the  wing  mean 
plane.  However,  a  modified  form  of  the  full  surface  boundary  condition  is 
used  in  place  of  the  small  disturbance  condition.  The  full  potential  boun¬ 
dary  condition  can  be  written  as 

♦*  "  (1  +  V  S  +  (V  dy  (15) 

In  the  current  design  method,  the  spanvise  slope  is  set  to  sero.  As 
in  the  small  disturbance  method,  $z  is  computed  using  Eq.  7  after  the 
scheme  has  Iterated  to  a  desired  level  of  convergence. 
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INVERSE  METHOD  VERIFICATION 


The  verification  teat  centered  on  validating  the  design  consistency 
and  accuracy  of  both  inverse  methods.  The  design  consistency  means  that  a 
wing  shape  generated  for  a  given  pressure  distribution  in  the  inverse  mode 
will  yield  the  same  pressure  distribution  when  run  in  a  purely  analysis 
mode.  A  test  of  both  consistency  and  accuracy  is  to  take  the  pressure 
distribution  for  a  known  wing  shape  as  the  target  pressure  distribution  for 
an  Inverse  design  and  then  compare  the  computed  shape  with  the  original 
wing.  Both  techniques  were  used  in  the  present  research. 

Small  Disturbance  Method 

Both  the  small  disturbance  and  the  full  potential  design  methods  were 
verified  using  an  untapered  NACA  0012  wing  with  an  aspect  ratio  of  6.96.  A 
planform  view  of  this  wing  configuration  Indicating  the  span  stations  used 
in  the  design  testing  is  shown  in  Figure  2.  Tests  were  made  for  both 
subcritical  and  supercritical  Mach  numbers.  All  the  tests  were  made  using 
a  72x17x30  grid. 

Figures  3  and  4  present  results  for  the  small  disturbance  design 
scheme  at  a  single  span  station  at  two  angles  of  attack.  These  figures 
compare  the  computed  pressures  and  target  pressures  for  two  subcritical 
tests.  In  Figure  3,  the  input  pressures  for  the  inverse  scheme  were  the 
same  as  those  obtained  from  analysis  for  a  2  degree  angle  of  attack.  As 
can  be  seen,  the  resulting  pressures  obtained  at  the  end  of  the  Inverse 
cycle  compare  well  with  the  analysis  pressures. 

Figure  4  presents  results  for  a  test  at  zero  angle  of  attack  using  a 
modified  upper  surface  pressure  distribution.  The  pressure  distribution 
obtained  at  the  end  of  the  Inverse  cycle  is  compared  with  the  target  pres- 
sure  distribution.  The  pressures  in  the  inverse  region  compare  quite  well. 
However,  the  pressures  on  the  lower  surface  and  in  the  nose  region  of  the 
upper  surface  are  slightly  changed.  These  changes  are  to  be  expected  since 
the  lift  and,  therefore,  the  circulation  of  the  wing  has  been  changed. 
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The  result*  obtained  for  the  emsll  dlsturbsnce  scheme  Indicated  that 
the  method  was  convergent  and  duplicable  for  a  given  pressure  distribu¬ 
tion.  Based  on  these  results,  the  method  wae  extended  to  the  full  poten¬ 
tial  scheae.  A  series  of  tests  were  run  in  which  the  inverse  scheae  was 
applied  at  three  consecutive  span  stations. 

Subcritical  Tests 

The  first  runs  with  the  full  potential  code  were  aade  to  test  the 
accuracy  of  the  code  by  uaing  the  pressure  distribution  for  a  known  airfoil 
shape  as  the  target  for  the  inverse  aode,  and  then  comparing  the  resulting 
shape  with  the  original  airfoil.  Figures  5  to  7  present  results  at  a 
single  span  station  for  a  subcritical  test  at  2  degrees  angle  of  attack. 
Figure  S  demonstrates  that  the  full  potential  aethod  will  accurately 
reproduce  the  desired  pressure  distribution.  The  accuracy  of  the  inverse 
scheme  and  the  shape  calculation  ia  shown  by  Figures  6  and  7.  A a  can  be 
seen,  both  the  KACA  0012  slopes  and  ordinates  are  accurately  computed. 

Figures  8  through  19  present  results  at  all  three  design  stations  for 
a  subcritical  test  at  zero  angle  of  attack  for  a  modified  upper  surface 
pressure  distribution.  The  target  pressure  distribution  was  obtained  by 
modifying  the  analysis  pressures  with  a  french  curve.  It  was  not  expected 
that  specifying  the  target  pressure  distribution  in  this  arbitrary  manner 
would  produce  realistic  airfoil  shapes.  However,  the  object  of  these  tests 
was  to  determine  if  the  specified  pressure  distribution  could  be  reproduced 
by  the  inverse  scheme. 

Figures  8  to  10  present  the  pressure  distribution  for  all  three  span 
stations.  As  previous  tests,  the  computed  pressures  agree  well  with  the 
specified  pressures  in  the  design  region  with  small  changes  evident  on  the 
lower  surface  and  the  nos*  region  of  the  upper  surface.  The  airfoil  elopes 
and  ordinates  computed  for  the  specified  design  pressures  are  given  in 
Figures  11  through  16.  Note  that  these  new  airfoils  all  have  "fishtails,” 
l.e.,  the  upper  and  lower  surfaces  cross,  which  is  physically  unrealistic. 
However,  these  results  illustrate  one  of  the  problems  encountered  in  devel¬ 
oping  an  inverse  scheme.  Because  the  purpose  of  these  tests  was  to  verify 
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Che  Inverse  approach  and  noc  co  design  a  wing,  no  attempt  was  made  to  con¬ 
trol  the  trailing  edge  thickness* 


Figures  17  through  19  illustrate  the  consistency  of  the  design  scheme. 
The  pressures  generated  by  analyzing  the  modified  wing  shape  are  compared 
with  the  pressures  obtained  from  the  inverse  code.  These  results  indicate 
that  the  modified  wing  shapes  will  yield  the  desired  pressures  when  the 
wing  is  analyzed. 

Supercritical  Tests 

The  next  series  of  tests  were  made  to  verify  the  Inverse  scheme  for 
supercritical  flow.  Following  the  same  procedures  used  for  subcritlcal 
flow,  the  code  was  first  tested  for  accuracy  by  using  a  pressure  distribu¬ 
tion  obtained  from  analysis  of  the  wing,  as  the  target  pressure  distribu¬ 
tion  for  the  inverse  scheme.  Figure  20  to  22  presents  results  for  the 
unswept  NACA  0012  wing  at  a  2  degree  angle  of  attack  and  a  Mach  number  of 
0.82.  As  in  the  subcritlcal  cases,  the  pressures  are  in  excellent 
agreement. 

The  next  step  in  the  supercritical  testing  was  to  modify  the  Mach  0.82 
two-degree  angle  of  attack  pressure  distribution  to  eliminate  the  shock  on 
the  upper  surface.  This  type  of  design  represents  a  typical  application  of 
an  inverse  scheme  in  a  wing  design.  Figures  23  to  25  compare  the  pressure 
distributions  obtained  by  the  inverse  scheme  with  the  desired  pressure 
distributions  at  the  three  design  stations.  These  results  further  confirm 
that  the  inverse  scheme  is  successful  at  supercritical  Mach  numbers. 

Figures  26  through  31  show  the  corresponding  airfoil  slopes  and  shapes 
obtained  from  the  Inverse  code.  Note  that  for  this  design  the  trailing 
edge  thicknesses  of  the  resulting  airfoils  are  physically  unrealistic. 
However,  Figure  32  illustrates  that  these  shapes  are  consistent  with  the 
specified  pressure  distribution.  The  results  at  the  other  two  design 
stations  compare  equally  well. 
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Swept  Wings 


A  final  series  of  tests  was  conducted  to  verify  the  inverse  method  for 
swept  wings.  The  supercritical  tests  were  repeated  using  the  basic  NACA 
0012  wing  swept  back  15  degrees.  The  Mach  number  for  the  tests  was 
increased  to  0.85.  The  results  presented  are  for  an  angle  of  attack  of 
zero  degrees.  Figures  33  to  35  compare  the  computed  and  specified 
pressured  distributions  for  a  test  in  which  the  upper  surface  pressures 
were  modified  to  eliminate  the  shock.  The  airfoil  slopes  and  shapes 
generated  by  the  modified  pressure  distribution  are  given  in  Figures  36 
through  38.  These  figures  show  that  the  accuracy  of  the  code  is  not 
affected  by  wing  sweep.  Figures  42  to  44  show  that  the  consistency  of  the 
design  scheme  is  retained. 

All  the  test  cases  described  were  converged  to  a  residual  of  0.0011, 
which  represents  a  reduction  from  the  initial  residual  of  about  three 
orders  of  magnitude.  For  most  cases,  this  represents  a  sufficient  level  of 
convergence.  This  convergence  criteria  required  an  average  of  about  83 
seconds  on  a  Cyber  203  for  inverse  runs  and  about  63  seconds  for  analysis 
runs. 


DEVELOPMENT  OF  THE  BASELINE  UNIFIED  DESIGN/ANALYSIS  CODE 

» 

The  baseline  unified  design/ analysis  code  was  developed  in  two  stages. 
In  the  first  stage,  a  three-dimensional  transonic  flow  analysis  code  was 
developed  that  served  as  the  foundation  for  the  development  of  the  baseline 
design/ analysis  code.  In  the  second  stage,  the  inverse  design  method 
developed  at  Texas  A&M  University  was  implemented. 

TRANSONIC  POTENTIAL  FLOW  ANALYSIS 

A  new  analysis  code  had  to  be  written  because  the  NASA  Langley  ZEBRA 
II  code  was  developed  as  a  pilot  code  to  verify  the  ZEBRA  algorithm  and 
was,  therefore,  not  suited  for  calculating  flow  about  arbitrary  wing  con¬ 
figurations.  Additionally,  the  ZEBRA  II  code  (and  consequently,  the  TAMU 
pilot  code)  was  written  in  explicit  vector  instructions  which  prohibits  its 
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use  on  computers  other  than  the  CYBER  203  and  205.  The  baseline  unified 
design/ analysis  code  was  written  in  standard  FORTRAN  to  make  it  transport¬ 
able  to  other  computers.  The  major  difference  between  the  baseline  code 
and  the  TAMU  code  is  the  use  of  streamwise  shearing  transformations  to 
align  the  computational  mesh  with  the  wing  planform. 

Grid  Generation 

Prior  to  the  start  of  this  research,  it  was  felt  that  the  fine  inner 
mean/ coarse  outer  mesh  grid  embedding  scheme  developed  by  Boppe^  would 
provide  the  optimum  mesh  system  for  the  design  code.  However,  after  work 
on  the  design  method  began,  it  was  decided  that  the  time  required  to  imple¬ 
ment  the  embedded  grid  system  would  delay  the  development  of  the  unified 
analysis/design  method.  Therefore,  both  analysis  and  inverse  solutions 
were  obtained  using  one  mesh  for  the  entire  computational  domain.  Initial 
tests  were  made  using  a  90x30x30  grid.  Because  the  convergence  rate  for 
this  grid  system  was  unacceptably  slow,  grid  sequencing  was  employed  to 
speed  up  convergence. 

The  computational  grid  for  the  baseline  analysis  code  was  formed  by 
first  computing  a  stretched  cartesian  grid  system  and  then  shearing  the 
grid  to  align  it  with  the  leading  and  trailing  edges  of  the  wing.  This 
procedure  has  been  used  with  great  success  in  the  small-disturbance  codes 
developed  by  Bailay  and  Ballhaua^  and  by  Boppe^.  With  the  sheared  grid 
system,  each  spanwise  plane  of  the  grid  contains  an  equal  number  of  points 
on  or  adjacent  to  the  wing  surface. 

To  ensure  compatibility  with  the  TAMU  pilot  code,  the  grid  stretching 
used  in  that  code  was  Implemented  in  the  Lockheed  code.  In  this  mesh 
system,  the  wing  surface  is  covered  with  an  evenly  spaced  grid  system.  The 
regions  in  front  of  the  leading  edge,  behind  the  trailing  edge,  and 
outboard  of  the  wing  tip  are  stretched  exponentially.  Geometric  stretching 
is  used  above  and  below  the  wing  mean  plane. 
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A  shearing  transformation  of  the  form 


X(x,y) 


x-xle(y) 

c(y) 


Y(y)  -  y 


Z(z)  -  z 


(16) 


(where  xle(y)  defines  the  leading  edge  of  the  wing  and  c(y)  is  the  local 
chord  distribution),  transforms  the  physical  grid  system  (x,y,z)  into  a 
computational  grid  (X,Y, Z)  aligned  with  the  wing.  Figure  45  shows  the  mesh 
system  generated  by  this  transformation. 

In  a  typical  inverse  or  analysis  solution,  a  sequence  of  three  grids 
is  used.  The  solution  starts  on  a  coarse  grid  that  has  25  chordwise 
points,  30  spanwlse  points,  and  8  points  normal  to  the  wing  mean  plane. 
The  solution  from  the  coarse  grid  is  interpolated  onto  a  medium  grid  that 
contains  50  chordwise  points,  30  spanwlse  points, and  16  normal  points.  The 
medium  grid  solution  is  interpolated  onto  a  fine  grid  that  has  90  chordwise 
points,  30  spanwlse  points,  and  30  points  normal  to  the  mean  plane.  Twenty 
of  the  30  spanwlse  stations  used  in  each  grid  were  placed  on  the  wing  sur¬ 
face.  The  number  of  chordwise  points  covering  the  local  chord  at  each  span 
station  varied  from  10  for  the  coarse  grid  to  25  for  the  medium  grid  and  50 
for  the  fine  grid. 


The  Full  Potential  Equation  in  General  Coordinates 

Following  Holst^  ®  a  general  coordinate  transformation  is  used  to 
transform  Eq.  (A-l)  to  the  computational  coordinate  system.  On  trans¬ 
formation,  Eq.  (A-l)  becomes 


p  -  [1  -  (U$x  +  V$Y  +  W*z)  y  '  1 1  (17b) 
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where  U,  V,  and  W  are  the  contravariant  components  of  velocity  in  the  com¬ 
putational  plane  and  J  is  the  Jacobian  of  the  transformation.  Details  of 
this  transformation  are  given  in  References  4  to  6.  Eq.  (17)  has  been 
nondimensionalized  by  the  critical  speed  of  sound  and  the  stagnation 
density. 

For  the  transformation  defined  by  Eq.  (16),  we  get 
U  -  (Xx2  +  Xy2)  *x  +  Xy 

V  -  Xy  +  ♦y 

(18) 

W  “  *Z 
J  “  X* 

This  transformation  retains  the  strong  conservation  form  of  the  original 
equation. 


Numerical  Solution  Algorithm 


The  finite  difference  analog  of  Eq.  (17)  can  be  written 

~  /£U\  +  $  (£l)  +  M2?) 

x  j)i^,J,k  Y 


(19) 


where  5X,  &  y  and  ^  are  first  order  backwards  differences.  In  order  to 
maintain  stability  in  regions  of  supersonic  flow,  the  density  has  been 
replaced  by  the  retarded  density  approximation  used  by  Holst. 


pi+4,j,k 


t(l-v)p]i+4>jfk+  vi+'i,j  ,k°i-4,j 


(20a) 


where 


M 

v  -  MIN  (1,  MAX  (1 - -  ,  0)] 

m2 

i.j.k 


(20b) 
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and  p  is  che  retarded  density  coefficient;  Mc  is  a  cutoff  Mach  number 
whose  value  is  usually  0.94<MCS1.0. 

In  the  current  code,  Eq.  (20)  is  evaluated  only  at  the  midsegment 
point  1+1/2, j,k.  The  values  at  i,j+l/2,k  and  i,j,k+l/2  are  obtained  by 
averages  of  the  surrounding  points.  Values  on  the  mean  plane  are  obtained 
by  two  point  extrapolation  from  above  and  below  the  mean  plane.  Averages 
of  central  difference  approximations  are  used  to  compute  the  values  of 
at  iil/2  and  at  J±l/2  in  Eq.  18. 

It  was  found  chat  the  convergence  of  the  analysis  method  was  improved 
by  splitting  the  full  potential  into  separate  perturbation  and  freestream 
components. 

4>  ■  G  +  X  q  cos  a  +  Zq  sin  a  (21) 

^oo  ao 


where  at  is  the  angle  of  attack.  Equation  (19)  is  then  solved  for  G. 
Boundary  Conditions 

As  in  the  NASA  Langley  ZEBRA  code,  surface  boundary  conditions  are 
introduced  by  replacing  W  at  KWNGB+1/2  and  KWNGT-1/2  in  Eq.  (19)  with  Eq. 
(15).  This  boundary  condition  is  computed  prior  to  the  start  of  each 
iteration  using  values  of  $  from  the  previous  iteration  and  then  held 
constant  during  the  ZEBRA  sweeps. 

For  lifting  cases,  the  section  circulation  r  is  computed  by  taking  the 
difference  in  the  potential  at  the  section  trailing  edge  linearly  extrapo¬ 
lated  from  the  points  above  and  below  the  airfoil.  The  Kutta  condition  is 
implemented  in  Eq.  (19)  by  replacing  4>kMNGb  with  4>KWHGB_r  on  the  line  above 
the  mean  plane  and  by  replacing  with  +r  on  the  llne  b*low  th< 

wing  plane  at  all  points  behind  Che  trailing  edge  of  the  wing.  The  circu¬ 
lation  is  computed  prior  to  the  start  of  each  iteration  and  then  slightly 
overrelaxed  to  improve  convergence. 
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Because  the  outer  boundaries  are  located  at  finite  lengths  from  the 
wing,  the  expressions  of  Klunker^  are  used  to  compute  the  change  in  the 
far-field  potential  due  to  lift.  The  potential  on  the  downstream  boundary 
is  updated  by  assumming  that  Gx  ■  0.  This  was  implemented  by  letting  - 

GNI-1 * 

The  symmetry  condition  at  the  wing  root  is  implemented  by  setting 
V  -  0  and  by  replacing  in  Eq.  (18)  with 

*Y  “  "  Xy  *X  (22) 

Verification  of  the  Analysis  Code 

A  standard  wing  used  to  evaluate  the  performance  of  an  analysis  code 
is  the  ONERA  M6  wing  described  in  Reference  12.  This  wing  has  a  leading 
edge  sweep  angle  of  30  degrees,  an  aspect  ratio  of  3.8,  and  a  taper  ratio 
of  0.562.  A  series  of  verification  tests  were  conducted  at  a  Mach  number 
of  0.839  and  an  angle  of  attack  of  3.0  degrees.  Figures  46  to  49  compare 
the  computed  pressure  distribution  at  four  span  stations  on  the  wing  with 
experimental  data  and  results  from  TWING  program  of  Holst.  As  can  be  seen, 
the  results  for  this  test  case  compare  reasonably  well  with  both  the  exper¬ 
imental  data  and  the  TWING  results.  The  discrepancies  at  the  wing  tip  and 
at  the  lower-surface  leading  edge  indicate  the  need  for  finer  grid  systems 
in  these  regions.  The  pressure  distribution  over  the  entire  wing  is  given 
in  Figure  50. 

In  a  typical  analysis  run,  the  coarse  and  medium  grids  are  Iterated 
until  the  initial  residual  drops  by  four  orders  of  magnitude  or  a  maximum 
number  of  iterations,  usually  400,  is  reached.  The  code  is  then  run  for 
100  iterations  on  the  fine  grid.  The  current  version  of  the  baseline  code 
takes  about  90  seconds  of  CPU  time  on  a  CRAY  IS  computer  and  about  8 
minutes  of  CPU  time  on  the  CYBER  203  to  run  the  ONERA  M6  test  case.  This 
disparity  in  run  times  points  out  the  inefficiency  of  the  automatic  vector- 
lzation  feature  of  the  CYBER  203  FORTRAN  compiler. 
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INVERSE  METHOD 


The  inverse  method  developed  at  Texas  A&M  University  was  implemented 
in  the  same  manner  as  was  done  in  the  ZEBRA  II  pilot  code.  However,  the 
non-dimensional  form  of  the  governing  equations  used  in  the  Lockheed  code 
leads  to  a  slightly  different  equation  for  the  pressure  coefficient 
boundary  condition.  The  derivation  of  this  equation  is  given  in  Appendix 
B. 


A  series  of  tests  were  initiated  to  determine  the  consistency  and 
accuracy  of  the  inverse  method  using  the  ONERA  M6  wing  as  the  base  geome¬ 
try.  First,  the  pressure  distribution  for  the  ONERA  M6  wing  generated  by 
analysis  for  the  Mach  0.84  case  was  used  as  the  target  pressure  distribu¬ 
tion  for  the  inverse  code.  The  resulting  pressure  distribution  and  wing 
shapes  were  then  compared  with  the  corresponding  data  for  the  base  shape. 
The  inverse  method  was  applied  at  five  consecutive  design  stations  along 
the  span.  Figures  51  to  53  compare  the  base  airfoil  sections  with  the 
sections  computed  by  the  inverse  scheme.  The  airfoil  shapes  are  recovered 
reasonably  well.  The  error  in  the  shape  calculation  is  felt  to  be  due  to 
neglecting  the  effects  of  the  spanwise  variation  in  slope  in  the  inverse 
scheme  for  a  tapered  wing.  In  addition,  the  target  pressures  were 
generated  using  the  full  potential  boundary  condition.  Techniques  to 
Include  the  spanwise  effects  in  the  inverse  scheme  will  be  investigated  in 
the  second  phase  of  this  contract. 

Next,  the  Inverse  scheme  was  tested  using  target  pressure  distribu¬ 
tions  designed  to  weaken  the  shock  at  each  of  the  five  design  stations. 
Figures  54  to  56  compare  the  computed  pressure  distributions  with  the 
target  pressure  distributions  at  three  of  the  design  stations.  The 
computed  and  the  target  pressures  are  in  good  agreement.  Figures  57  to  59 
compare  the  base  airfoils  with  the  computed  airfoils  by  the  inverse  method. 
As  in  the  results  obtained  by  the  TAMU  pilot  code,  the  airfoils  produced 
are  physically  unrealistic.  Figures  60  to  62  Indicate,  however,  that  the 
computed  wing  shape  is  consistent  with  the  specified  pressure  distribution. 
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In  a  typical  inverse  run,  the  coarse  grid  is  converged  for  the  base 
wing  geometry.  The  coarse  grid  results  are  interpolated  onto  the  medium 
grid  and  run  for  SO  iterations  before  the  inverse  design  is  initiated.  The 
medium  grid  is  then  run  in  the  inverse  mode  until  convergence.  The  medium 
grid  results  are  interpolated  onto  the  fine  grid  which  is  then  run  for  the 
100  iterations  in  the  inverse  mode.  Wing  shapes  are  computed  at  the  end  of 
the  fine  grid  iterations.  Running  the  code  in  the  inverse  mode  requires 
only  about  100  seconds  on  the  CRAY  IS. 

CONTROL  OF  TRAILING  EDGE  THICKNESS 

The  results  from  both  the  baseline  and  pilot  inverse  design  codes 
indicate  the  need  for  a  technique  to  control  trailing  edge  thickness.  A 
review  of  existing  wing  and  airfoil  Inverse  design  methods  revealed  two 
different  procedures  for  enforcing  trailing  edge  closure  that  can  be  used 
in  the  present  design  method. 

In  the  first  approach,  the  nose  region  of  the  starting  airfoil  or  wing 
is  modified  to  increase  or  decrease  the  leading  edge  radius.  After  a  few 
tries,  a  nose  shape  can  usually  be  found  that  will  provide  the  desired 
trailing  edge  thickness.  This  approach  was  used  by  Carlson13-14  in  his 
two-dimensional  direct/ inverse  design  method.  The  major  drawback  of  this 
approach  is  that  it  is  a  trail  and  correction  procedure  that  relies  heavily 
on  the  expertise  of  the  designer.  Shankar13-1^  has  suggested  that  this 
procedure  can  be  automated  by  specifying  the  nose  shape  by  y  •  a0xn,  where 
n  and  a  are  free  parameters  that  are  adjusted  by  a  numerical  optimization 
procedure  to  satisfy  a  specified  trailing  edge  thickness  constraint. 

A  second  technique  to  enforce  trailing  edge  closure  has  been  used  by 
Shankar1^  in  a  small-disturbance  design  code.  In  this  technique,  a 
functional  relationship  between  the  trailing  edge  thickness  and  the 
velocity  potential  at  the  leading  edge  is  assumed  at  each  spanwlse  design 
station.  A  perturbation  to  the  leading  edge  velocity  that  will  drive  the 
trailing  edge  thickness  to  a  desired  value  can  then  be  computed  at  each 
design  station. 
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At  the  start  of  this  contract  it  was  felt  that  perturbing  the  nose 
shape  under  control  of  a  numerical  optimization  scheme  would  prove  to  be 
the  most  effective  approach  for  enforcing  trailing  edge  closure. 
Therefore,  the  baseline  unified  design/ analysis  code  was  modified  to  use 
the  optimization  techniques  developed  by  Vanderplaats,  Hicks,  and  Murman*7 
to  perform  the  nose  shape  modifications  required  to  control  trailing  edge 
thickness.  However,  the  time  constraints  on  the  first  phase,  of  the 
contract  prevented  the  successful  implementation  of  optimization.  The  best 
way  to  control  the  trailing  thickness  will  be  addressed  in  the  second  phase 
of  the  contract. 


RESULTS  AND  DISCUSSION 

The  results  presented  in  the  previous  sections  of  this  report 
illustrate  that  the  present  design  method  can  be  used  effectively  for  both 
subcrltical  and  supercritical  design  cases.  The  supercritical  results  show 
that  the  method  can  perform  one  of  the  more  important  functions  of  a 
transonic  wing  design  code— the  elimination  or  weakening  of  strong  shock 
waves  at  supercritical  Mach  numbers.  This  can  be  seen  by  comparing  the 
section  lift,  drag,  and  moment  coefficients  shown  in  Fig.  48  for  the  67Z 
span  station  of  the  base  ONERA  wing  with  the  values  given  in  Fig.  62  for 
modified  wing.  The  lift  coefficient  changes  from  .21  for  the  base  wing  to 
.020  for  the  modified  wing.  The  drag  coefficient  changes  from  -0.0081  to 
-0.001  and  the  moment  coefficient  changes  from  -0.0960  to  -0.0947.  The 
negative  values  of  drag  illustrate  the  Inaccuracy  of  the  pressure  drag 
calculation  that  is  common  to  most  transonic  analysis  codes.  However,  the 
values  obtained  illustrate  the  effect  that  weakening  the  shock  has  on  the 
section  drag. 


CONCLUSIONS 

An  Inverse  design  method  for  wings  in  transonic  flow  that  is  particu¬ 
larly  suited  for  use  on  a  vector  computer  has  been  developed.  The 
technique  has  been  verified  for  accuracy  and  consistency  in  both  a 
developmental  pilot  code  and  a  baseline  unified/design  code  that  will  serve 
as  the  basis  for  future  code  development. 
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Figure  5.  Correlation  of  full  potential  inverse  and  analysis 
target  pressures,  Mach  ■  0.4, n  “  0.3125 


Figure  7 .  Correlation  of  original  airfoil  ordinates  with  ordinates 
from  inverse  using  analysis  target  pressures 


Figure  8.  Comparison  of  inverse  pressures  and  modified  target 
pressures,  Mach  ■  0.4,  n  ■  0.1875 


30 


Figure  10.  Comparison  of  inverse  pressures  and  modified  target 
pressures,  Mach  *  0.4, n  ■  0.4375 
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Figure  13.  Airfoil  alopea  obtained  from  modified  target  preaaurea 
Mach  -  0.4.  n  -  0.3125 
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Figure  16.  Airfoil  ordinates  obtained  from  modified  target  pressurea 
Mach  -  0.4,  n  -  0.4375 


38 


-0.50 


-0.45 


-0.40 


-003  4. 


-0.30 


-0.25 


-0.20-1 


Straight  ting 
A2»6.96 
M  *0.4 
a  *0.0 

Wing  Station  3 


Lower  Surface 


Uocer  Surface 


-0.15 ■ 


-0.10- 


-0.05- 


Figure  18.  Correlation  of  inverse  pressures  and  pressures  from 
analysis  of  modified  wing,  Mach  ■  0.4,  n  *  0.3125 


40 


Straight  Wing 
AS*6.56 

M  »C.4 

3  =  0.0 

Wing  Station  5 

\ 

\ 

Uooer  Surface 

over  Surface 

A 

V 

\ 

- 1 - 1— 

|  0.1  0.2 

)  1 - 1 - 1 - -  >  \  \  - - - — 

0.3  0 .4  0.5 _ 0.6  0.7  o.av\o.9  1.0 

-Topj  Prsaaurea  Proa  inverse 

0 

gsass’AiS&r1 37  \v 

.Sottoai  Preaaurea  Proa  Inverse  Hi 

A 

Pressures  Cn  Sot ten  When  \\ 

Analysing  Modified  Airfoil  \\^ 

-* 

k 

Figure  19.  Correlation  of  inverse  pressures  and  pressures  from 
analysis  of  modified  wing,  Mach  “  0.4,  n  “  0.4375 


Figure  20.  Correlation  of  inverae  and  analysis  target  pressures 
Mach  -  0.82,  *  -  2.0,  n  -  0.1875 
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Figure  21.  Correlation  of  inverse  and  analysis  target  pressures, 
Mach  ■  0,S?.,a-  2.0,  n  -  0.3125 


Figure  22,  Correlation  of  inverae  end  analyaia  target  preaaurea, 
Mach  -  0.82,  2.0,  n  -  0.4375 
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23.  Correlation  of  inverse  and  modified  target  pressures 
Mach  ■  0.82,  a  -  2.0,  n  ■  0.1875 


Figure  24.  Correlation  of  inverse  and  modified  target  pressures, 
Mach  -  0.82,  3“  2.0,  n  “  0.3125 
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Figure  25.  Correlation  of  inverse  and  aodified  target  pressures, 
Mach  -  0.82,  a  -  2.0,  n  •  0.4375 
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Figure  26.  Airfoil  >lope>  from  modified  target  preaiuree, 
Mach  -  0.82,  a*  2.0,  n  •  0.1875 
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Figure  29.  Airfoil  ordinates  from  modified  target  pressures, 
Mach  -  0.82,  a  -  2.0,  n  -  0.3125 
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Figure  32.  Comparison  of  inverae  preaaures  with  pressures  from 
analysis  of  modified  wing,  n  ■  0.1875 
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Figure  33.  Correlation  of  modified  target  presaurea  and  inverae 
preaaurea  for  a  awept  wing,  Mach  ■  0.85,  *  *  2.0, 
n  -  0.1875 
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Figure  35.  Correlation  of  modified  target  pressures  and  inverse 
pressures  for  a  swept  wing,  Mach  *  0.85,  a  *  2.0, 
n  -  0.4375 
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Figure  38.  Comparison  of  inverae  pressures  with  pressures  from 

analysis  of  modified  swept  wing,  Mach  "  0.85,  a  -  2.0, 
n  -  0.4375 


Figure  39.  Modified  avept  wing  eirfoil  slopes,  n  ■  0.1875 
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Figure  41.  Modified  swept  wing  airfoil  slopes,  n  ■  0.312$ 


Figure  42.  Modified  avept  wing  airfoil  ordinatea,  n  ■  0.3125 
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Figure  45.  Planfora  view  of  aheared  cartesian  grid 
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Figure  46.  Correlation  of  pressures  from  analysis  of  0NEBA  M6  wing 
with  experiment  and  TWING,  n  -  0.20 
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Figure  47.  Correlation  of  pressures  from  an* lysis  of  ONERA  M6  wing 
with  experiment  end  TWING,  n  -  0.4103 
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Figure  48.  Correlation  of  preaaurea  from  analysis  of  ONERA  M6  ving 
with  experiment  and  TWXNG,  n  -  0.6667 


Comparison  of  ONERA  M6  ordinates  with  ordinates  from 
inverse  using  analysis  target  pressures,  n  ”  0.4615 
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Figure  54.  Comparison  of  inverse  pressures  with  modified  target 
pressures,  n  -  0.4615 
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Figura  55.  Comparison  of  invars#  pressures  with  modified  target 
pressures,  n  -  0.5641 
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Figure  56.  Conpariaon  of  inverae  preeeurea  with  modified  target 
preaaurea,  n  -  0.6667 
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Comparison  of  ONERA  M6  airfoil  ordinates  and  ordinates 
from  inverse  using  modified  target  pressures 9  n  ■  0#56« 


Comparison  of  ONERA  M6  airfoil  ordinates  and  ordinates 
from  inverse  using  modified  target  pressures,  n  ■  0.4615 


Comparison  of  ONERA  H6  airfoil  ordinatea  and  ordinates 
from  inverse  using  modified  target  pressures,  n  ■  O.661 


MRCH- 
CL  - 


0.8395  RLPHR-  3.0600  TWIST-  0.0000 

0.2041  CM  --0.0947  CD  --0.0010 


Figure  62.  Correlation  of  inverse  pressures  and  pressures  from 
analysis  of  modified  wing,  n  “  0.6667 


APPENDIX  A 
ZEBRA  II  ALGORITHM 


The  ZEBRA  II  scheme  of  South  et  al.^-^  solves  the  conservative  form  of 
the  full  potential  equation 


where 


(pu)  +  (pv)  +  (pw)  *  0 
x  y  z 


u  *  <)> 


v  “  4> 


(A-  1  ) 


w  “  $ 
p  =*  (M2  a2) 


z  1 

Y-l 


(  A-2  ) 


a2  -  —  +  I^i(1  _  q2) 
M2  2 


J2  «  u2  +  V2  +  w2 


Equation  (A-l)  is  replaced  by  its  finite  difference  analog  on  an 
evenly  spaced  grid 


5X  <P*X>  .  .+  dy  (0*y>  ^  +  6Z  *  0  (A-3) 

i+H.j.k  7  7  i,j-^S,k  2  z  i ,  j  ,  k+4 


where  5^,  6  ,  and  i>z  are  first-order  backwards  differences.  For  example, 


5x  <0*x> 


where 


i+h ,  j  ,k 


(.oil  )  -  (oi>  ) 

1+4. ik  x  i-^.ik 
Ax 


(A-4) 


i+4,jk  Ax 


(A-5) 


$ 


x  „  0i,jk  ~  *i-l,1k 

i-*a,Jk  Ax 


The  density  o  is  replaced  by  the  upwinded  artificial  compressibility 
value  o  to  stabilize  calculations  In  regions  of  supersonic  flow.  The  value 
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of  p  is  given  by 


p  -  p  -  u  rJL  Ax  _3p_  _v_  Ay  3p,  (A_6) 

where  3o/3x  and  3c/py  are  upwind  differenced  at  supersonic  points,  Ax  is 
the  chordwlse  grid  spacing  and  Ay  is  the  spanwise  grid  spacing.  The 
upwind  switching  function  u  is  given  by 

a2 

V  -  MAX  (0,1 - )  (  A-  7  ) 

q2 

In  the  present  ZEBRA  II  code  v/q^  and  3p/3y  are  assumed  to  be  negligible 
and  u/q^,  is  assumed  to  be  approximately  one.  Equation  (A-6)  then  becomes 

3p  ,  „ . 

p  ”  p  -  u  r~  (A-8) 


The  ZEBRA  II  algorithm  solves  Eq.  (A-3)  using  an  iterative  scheme  that 
mimics  point  Successive  Overrelation  (SOR).  For  3-D  calculations,  the 
ZEBRA  II  algorithm  marches  in  the  streamwise  (I)  direction  solving  one 
spanwise  plane  at  a  time.  In  each  plane,  points  J+K  odd  are  denoted  black. 
and  points  J+K  even  are  termed  white.  Each  plane  is  solved  by  a  two-pass 
sweep  in  which  new  black  values  are  obtained  first,  followed  by  the  white 
points.  In  this  way,  convergence  is  accelerated  because  calculations  at 
the  white  points  will  use  updated  quantities  at  the  black  points. 

By  replacing  0  in  Eq.  (A-3)  with 

-  N+l  ,  M  (A_9) 


N+l  N 

where  AO  ”  i  -0  is  the  correction,  ^  is  an  acceleration  parameter,  and  N 
is  the  iteration  number.  The  solution  at  each  grid  point  is  given  by 


Ap 


l.j.k  B 


Ax2  R 


SAP 


ave 


i-l.j.k 


(  A-  1  0 ) 


8b 
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F/G  20/4 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS  196 a  A 


where 


B  -  -  +  -  (t*)2  +  —  (|f)2  +  6 

uj  Ay  u>  Az 
x  y  J  z 


(A-il) 


The  SAqj^  ^  term  In  Eq.  (A-10)  comes  from  the  Inclusion  of  a  term  to 

add  explicit  temporal  damping  to  the  algorithm;  B  is  the  damping  coeffic¬ 
ient*  It  should  be  noted  that  the  ZEBRA  II  algorithm  possesses  some  natu¬ 
ral  temporal  damping  since  the  points  in  the  1+1  plane  are  not  as  updated 
as  the  points  in  the  1-1  plane.  In  addition,  the  black  points  in  each  I 
plane  are  not  as  current  as  the  white  points. 


The  residual,  Ax2R.  ,  .  in  Eq.  (A-10)  can  be  written 

/  n  m  \  /  N  N+l  \ 

i.j.u  ■  •  ♦lj.k)"  + 


Ax2  R 


v  •  N  for  J+K  odd 
v  •  N+l  for  J+K  even 


(A-12) 
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APPENDIX  B 

DERIVATION  OP  PRESSURE  COEFFICIENT  BOUNDARY  CONDITION 
FOR  BASELINE  UNIFIED  DESIGN/ ANALYSIS  CODE 


Normalizing  the  full  potential  equation  by  the  critical  speed  of  sound 
and  the  stagnation  density  leads  to  a  slightly  different  equation  for  the 
pressure  coefficient  than  that  uaed  in  the  TAMU  pilot  code.  The  coef¬ 
ficient  of  pressure  can  be  written 


cp  ”  ~~  jr  -  U 

ym2  p8  p“ 


(B-l) 


,  where  P  is  the  local  pressure,  Ps  is  the  stagnation  pressure,  and  Poo  is 
the  free  stream  pressure.  Thus, 


-  [i  +  Jt!  m£]  Y~l 

Poo  2  “ 


(B-2) 


_P_ 

P 


(B-3) 


11  ‘  wT 


♦*+♦?!  Y'‘ 


(B-4  ) 


NOW  let  Pg / ?od 


and 


+ 


$2  $2 
*2  (1  +  _JL  +  _2) 

*1 


(B-5) 


Substituting  into  Eq.  (B-l)  yields 


cp  +  l  -  [l  -  *x  a  +  ^  +  — )!  Y"1  (B“6) 


2  .  Y-l 


2  $2 

X  X 
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Solving  Eq.  (B-6)  yields 


The  procedure  outlined  in  the  description  of  the  development  of  the 
three-dimensional  inverse  pilot  code  is  used  to  extract  a  value  of  <t>  as  a 
function  of  Cp  from  Eq.  (B-7  )  for  use  as  a  Dirichlet  boundary  condition  in 
Che  inverse  design  scheme. 
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